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Abstract 



This paper discusses a general formulation of the material point method in the context of additive 
decomposition rate-independent plasticity. The process of generating the weak form shows that volume 
integration over deforming particles can be major source of error in the MPM method. Several useful 
Q-i identities and other results are derived in the appendix. 

I" 

1 Governing Equations 

o 

^ The equations that govern the motion of an elastic-plastic solid include the balance laws for mass, 

momentum, and energy. Kinematic equations and constitutive relations are needed to complete the system 
of equations. Physical restrictions on the form of the constitutive relations are imposed by an entropy 

^h * inequality that expresses the second law of thermodynamics in mathematical form. 

The balance laws express the idea that the rate of change of a quantity (mass, momentum, energy) in a 
volume must arise from three causes: 

J> 1. the physical quantity itself flows through the surface that bounds the volume, 

2. there is a source of the physical quantity on the surface of the volume, or/and, 

l> 

3. there is a source of the physical quantity inside the volume. 

Let 17 be the body (an open subset of Euclidean space) and let <9$7 be its surface (the boundary of Q). 
Let the motion of material points in the body be described by the map 

> x = ¥ >(X) = x(X) 

X 

where X is the position of a point in the initial configuration and x is the location of the same point in the 
deformed configuration. The deformation gradient (F) is given by 

r, ^X „ 
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1.1 Balance Laws 



Let /(x, t) be a physical quantity that is flowing through the body. Let g(x, t) be sources on the surface of 
the body and let h(x., t) be sources inside the body. Let n(x, t) be the outward unit normal to the surface 
dfl. Let v(x, t) be the velocity of the physical particles that carry the physical quantity that is flowing. 
Also, let the speed at which the bounding surface dU is moving be u n (in the direction n). 
Then, balance laws can be expressed in the general form ([ 1]) 



d 
dl 



/(x,t) dV 



/ /(x,t)K(x,i)-v(x,i)-n(x,i)]dA+ / 5 (x,t)dA+ / h(x,t)dV 
Jan Jan Jo. 



Note that the functions /(x, t), ^(x, t), and /i(x, t) can be scalar valued, vector valued, or tensor valued - 
depending on the physical quantity that the balance equation deals with. 

It can be shown that the balance laws of mass, momentum, and energy can be written as (see Appendix): 

p + p V • v = Balance of Mass 

pv — V-cr — pb = Balance of Linear Momentum 

<x = <j Balance of Angular Momentum 

p e — cr : ( Vv) + V • q — p s = Balance of Energy. 

In the above equations p(x, t) is the mass density (current), p is the material time derivative of p, v(x, t) is 
the particle velocity, v is the material time derivative of v, tr(x, t) is the Cauchy stress tensor, b(x, t) is the 
body force density, e(x, t) is the internal energy per unit mass, e is the material time derivative of e, q(x, t) 
is the heat flux vector, and s(x, t) is an energy source per unit mass. 
With respect to the reference configuration, the balance laws can be written as[j] 

p det(F) — po = Balance of Mass 

po x — Vo • P T — po b = Balance of Linear Momentum 

F ■ P = P T ■ F T Balance of Angular Momentum 

rp 

po e — P : F + Vo • q — po s = Balance of Energy. 

In the above, P is the first Piola-Kirchhoff stress tensor, and po is the mass density in the reference 
configuration. The first Piola-Kirchhoff stress tensor is related to the Cauchy stress tensor by 

P = det(F) F 1 a . 

The gradient and divergence operators are defined such that 

Vv= ^ J^-ei (g) ej = Vijei (g) a,- ; V • v = ^ ^ = u M ; V • S = ^ = a ij:j ej . 

i,j=l ^ i=l 1 i,j=l ^ 



We can alternatively define the nominal stress tensor N which is the transpose of the first Piola-Kirchhoff stress tensor such 

that 

N ~ P T = det(F) (F _1 • <t) t = det(F) a ■ F~ T . 

Then the balance laws become 

p det(.F) — po = Balance of Mass 

po x — Vo ■ N — po b = Balance of Linear Momentum 

F ■ N T = N ■ F T Balance of Angular Momentum 

po e — N : F + Vo ■ q — po s = Balance of Energy. 
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where v is a vector field, BS is a second-order tensor field, and are the components of an orthonormal 
basis in the current configuration. Also, 



3 r\ 3 q 3 

where v is a vector field, i?5 is a second-order tensor field, and Ej are the components of an orthonormal 
basis in the reference configuration. 
The contraction operation is defined as 

3 

A .B = Ay j = Aj Bij . 
1.2 The Clausius-Duhem Inequality 

The Clausius-Duhem inequality can be used to express the second law of thermodynamics for 
elastic-plastic materials. This inequality is a statement concerning the irreversibility of natural processes, 
especially when energy dissipation is involved. 

Just like in the balance laws in the previous section, we assume that there is a flux of a quantity, a source of 
the quantity, and an internal density of the quantity per unit mass. The quantity of interest in this case is the 
entropy. Thus, we assume that there is an entropy flux, an entropy source, and an internal entropy density 
per unit mass (r/) in the region of interest. 

Let 17 be such a region and let dfl be its boundary. Then the second law of thermodynamics states that the 
rate of increase of rj in this region is greater than or equal to the sum of that supplied to $7 (as a flux or from 
internal sources) and the change of the internal entropy density due to material flowing in and out of the 
region. 

Let dU move with a velocity u n and let particles inside Q have velocities v. Let n be the unit outward 
normal to the surface dQ. Let p be the density of matter in the region, q be the entropy flux at the surface, 
and r be the entropy source per unit mass. Then (see HI [2j ), the entropy inequality may be written as 



d 
dt 



/ pr\ dV I > / prj (u n — v • n) dA + / q dA + pr dV 

Jn ) Jan Jan Jn 



The scalar entropy flux can be related to the vector flux at the surface by the relation q = — ^(x) • n. Under 
the assumption of incrementally isothermal conditions (see [3 ] for a detailed discussion of the assumptions 
involved), we have 

it \ q ( X ) S 
W) = r = - 

where q is the heat flux vector, s is a energy source per unit mass, and T is the absolute temperature of a 
material point at x at time t. 

We then have the Clausius-Duhem inequality in integral form: 



[pvM)>f Pv(u n -v-n)dA- f 5-^dA+ / ^JdV. 
Jn / Jan Jan 1 Jn 1 

We can show that (see Appendix) the entropy inequality may be written in differential form as 



dt 



prj > —V 
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In terms of the Cauchy stress and the internal energy, the Clausius-Duhem inequality may be written as (see 
Appendix) 

q-VT 

p (e - T 17) - <t : Vv < — . 

1.3 Constitutive Relations 

A set of constitutive equations is required to close to system of balance laws. For large deformation 
plasticity, we have to define appropriate kinematic quantities and stress measures so that constitutive 
relations between them may have a physical meaning. 

1.3.1 Thermoelastic Relations 

Let the fundamental kinematic quantity be the deformation gradient (F) which is given by 

^=^ = V x; detF>0. 

A thermoelastic material is one in which the internal energy (e) is a function only of F and the specific 
entropy (77), that i^] 

e = e{F,r]) . 

For a thermoelastic material, we can show that the entropy inequality can be written as (see Appendix) 

fde \ ( de „ T \ ■ q-VT 

P {d-,- T ) r]+ { P dF-- (T - F )-- F + — 

At this stage, we make the following constitutive assumptions: 

1. Like the internal energy, we assume that er and T are also functions only of F and 77, i.e., 

* = *(F,7 1 ); T = T(F,r 1 ). 

2. The heat flux q satisfies the thermal conductivity inequality and if q is independent of rj and F, we 
have 

q-VT<0 =>- -(k-VT)-VT<0 k>0 
i.e., the thermal conductivity k is positive semidefrnite. 



2 If a thermoelastic body is subjected to a rigid body rotation Q, then its internal energy should not change. After a rotation, the 
new deformation gradient (F) is given by 

F = Q F . 

Since the internal energy does not change, we must have 

e = e(F,77) = e(F,77). 

Now, from the polar decomposition theorem, F = R ■ U where R is the orthogonal rotation tensor (i.e., R ■ R T — R T ■ R — 1) 
and U is the symmetric right stretch tensor. Therefore, 

e(Q-R-U,n)=e(F, v ). 

Now, we can choose any rotation Q. In particular, if we choose Q — R T , we have 

e(R T ■ R ■ U, rj) = e(l ■ U,f]) = e(U,n) . 

Therefore, 

e(U, V ) = e(F, V ) . 

This means that the internal energy depends only on the stretch U and not on the orientation of the body. 
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Therefore, the entropy inequality may be written as 

>(M- T ) T ) - F< -°- 

Since rj and F are arbitrary, the entropy inequality will be satisfied if and only if 

de m m de , de „ T de „ T 

Therefore, 

Given the above relations, the energy equation may expressed in terms of the specific entropy as (see 
Appendix) 

pT rj = — V • q + p s . 
1.3.2 Alternative strain and stress measures 

The internal energy depends on F only through the stretch U. A strain measure that reflects this fact and 
also vanishes in the reference configuration is the Green strain 

E= 1 -(F T -F-1)= 1 -(U 2 -1). 



Recall that the Cauchy stress is given by 



de tpT 

a = P dF- F ■ 



We can show that the Cauchy stress can be expressed in terms of the Green strain as (see Appendix) 

„ de rp 
<r = P F.-.F-. 

Recall that the nominal stress tensor is defined as 

N = det F (a ■ F~ T ) . 
From the conservation of mass, we have po = P det F. Hence, 

N= P -^a.F- T . 

P 

The nominal stress is unsymmetric. We can define a symmetric stress measure with respect to the reference 
configuration call the second Piola-Kirchhoff stress tensor (S): 

S :=F- 1 -N = P F~ T = — F- 1 ■ a ■ F~ T . 

P 

In terms of the derivatives of the internal energy, we have 

*=f*-'-("-£^)-^=»£ 



and 



Po ( „ de T \ T ^ de 



That is, 

de de 
S = Po dE and N = PoF -dE 
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1.3.3 Stress Power 

The stress power per unit volume is given by a : Vv. In terms of the stress measures in the reference 
configuration, we have 

-:Vv=( /5j P.||.F^:(F.F- 1 ). 
Using the identity A : (B ■ C) = (A ■ C T ) : B, we have 



<r : Vv 



: F = p I F ■ ^ ) :F=—N:F 
dE) p 



We can alternatively express the stress power in terms of S and E. Taking the material time derivative of E 
we have 

E= ^(F T -F + F T -F) . 

Therefore, 

S : E = -[S : (F T ■ F) + S : (F T ■ F)} . 

Using the identities A: (B ■ C) = (A • C T ) : B = (B T • A) : C and A : B = A T : B T and using the 
symmetry of S, we have 

S : E = -[(S ■ F T ) : F T + (F ■ S) : F] = -[(F ■ S T ) : F + (F ■ S) : F] = (F ■ S) : F . 

Now, S = F^ 1 ■ N. Therefore, S : E = N : F. Hence, the stress power can be expressed as 

a :Vv = AT : F = S : E . 
If we split the velocity gradient into symmetric and skew parts using 

Vv = I = d + w 

where d is the rate of deformation tensor and w is the spin tensor, we have 

a : Vv = a : d + a : w = tr (er T ■ d) + tr (<r T ■ w) = tr (<r ■ d) + tr (<r ■ w) . 

Since cr is symmetric and w is skew, we have tr (cr ■ w) = 0. Therefore, cr : Vv = tr (er • d). Hence, we 
may also express the stress power as 



tr {tr ■ d) = tr (n t ■ = tr (s ■ E^j . 



1.3.4 Helmholtz and Gibbs free energy 

Recall that 

Therefore, 
Also recall that 



de 

S = p0 dE 



— = -S 
dE po 



de 



Now, the internal energy e = e(E, rj) is a function only of the Green strain and the specific entropy. Let us 
assume, that the above relations can be uniquely inverted locally at a material point so that we have 

E = E{S,T) and r] = fj{S,T). 

Then the specific internal energy, the specific entropy, and the stress can also be expressed as functions of S 
and T, or E and T, i.e., 

e = e(E,r)) = e(S,T) = e(E,T); rj = fj(S,T) = rj(E,T) ; and S = S(E,T) 
We can show that (see Appendix) 

4(e-T rj) = -f f]+ — S : E or ^ = -t rj + — S : E . 
at po at po 



and 



ie-Ti/-- S:E) = -fr,-— S : E or — =Tri+ — S: E. 
at po po at po 



We define the Helmholtz free energy as 

ip = ijj(E,T) :=e-Tr}. 

We define the Gibbs free energy as 

g = g(S,T) :=-e + Tr, + -S:E. 

Po 

The functions %[)(E, T) and g(S, T) are unique. Using these definitions it can be showed that (see 
Appendix) 



and 



dS df) dE dfj 

df = - po dE md df = p0 dS 



1.3.5 Specific Heats 

The specific heat at constant strain (or constant volume) is defined as 

_ de(E,T) 

The specific heat at constant stress (or constant pressure) is defined as 

de(S,T) 



We can show that (see Appendix) 



c P - dT 



v ~ dT~ or 2 
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and 

dfj 1 dE _ d 2 g d 2 g 



C v = T-^ + —S: — = T—^ + S: 



p dT po ' dT dT 2 ' dSdT ' 
Also the equation for the balance of energy can be expressed in terms of the specific heats as (see Appendix) 

pC v t = V ■(k-VT) + P s + — T0s :E 

Po 

p \C p -—S:a E ) f = V • (k ■ VT) + p s - — T a E : S 
\ Po J Po 

where 

PS '■= t^; and cxe '■= t^t • 
dT dT 

The quantity (3s is called the coefficient of thermal stress and the quantity ole is called the coefficient of 
thermal expansion. 

The difference between C p and C v can be expressed as 

C C - 1 is T di ]- d£! 

However, it is more common to express the above relation in terms of the elastic modulus tensor as (see 
Appendix for proof) 

1 T 

C p - C v = — S : cue H a E :C: a E 

Po Po 

where the fourth-order tensor of elastic moduli is defined as 

OS dH 
C := — - = po — = — - . 
dE dEdE 

For isotropic materials with a constant coefficient of thermal expansion that follow the St. Venant-Kirchhoff 
material model, we can show that (see Appendix) 

C p -C v = — [a tr (S) + 9 a 2 K T] . 
Po 

1.3.6 Kinematic Relations 

Additive split of the rate of deformation tensor: The rate of deformation tensor (d) is given by 

d= I[ V v + (Vv) T ] 

where v is the velocity. We assume that the rate of deformation can be additively decomposed into an 
elastic part, a plastic part, and a thermal part: 

d = d e + d p + d th . 

The thermal part of the rate of deformation is computed separately and subtracted from d to give 

d := d - d th = e e + dP . 
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For simplicity, we use the symbol d instead of d in the following development. 
We split the rate of deformation tensors into volumetric and deviatoric parts: 

d=\ei+rj; d e = \tx{d e ) 1 + V e ; dP = \ tr (tP) 1 + rf 

where 6 = tr (d), the deviatoric part of the rate of deformation tensor is r) = d — g tr (d) 1 , and 1 is the 
second-order identity tensor. 



1.3.7 Stress 

To deal with nearly incompressible behavior, we introduce an isomorphic split of the Cauchy stress tensor 
into a volumetric and a deviatoric part of the form (0): 

a = a m l+a s s a = p 1 + s (1) 

where 

1 r 1 1 

a m = — tr(cr) = V3p ; a s = \\s\\ ; p = - tr(<r) ; s = a - - tr (a) 1 

1 ^ s , 

1 = 71 — s = t: — — : A = V A : A . 
||s|| 

Note that the following identities hold: 

a m = a : 1 ; c s = er:s; 1:1=1; s : s = 1 ; 1 : s = ; 1 : s = ; s:s = 0. 

The rate form of equation ([T} is 

& = a m 1 + & s s + a s s & = p 1 + s . (2) 

In the following development, we enforce frame indifference while evaluating the constitutive relation by 
assuming that the stresses and rate of deformation have been rotated to the reference configuration using 

s = R T s R ; r} = R T r] R 

where the rotation tensor R is obtained from a polar decomposition of the deformation gradient. 



1.3.8 Elastic Relations: 

The constitutive model is assumed to consist of two parts. The first is a hypoelastic model for the deviatoric 
part of the stress and the second is a equation of state for the mean stress. 
The deviatoric stress is given by the rate equation 

s = C(a m ,T) :rf ^=> s = t(p,T):vf 

where s is the deviatoric stress rate, C is a pressure and temperature dependent fourth order elastic modulus 
tensor, and rf is the deviatoric part of d e (the elastic part of the rate of deformation tensor). 
Let us assume that the deviatoric elastic response of the material is isotropic, i.e. , [^] 

s = 2/i( C j m ,r)77 e ^ s = 2 Jl(p, T) rf . (3) 

3 Note that since the shear modulus depends on <r m , the rate of s should have a contribution that contains <r m . We ignore this 
contribution in the subsequent development. 
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The mean stress (<r m ) is calculated using a Mie-Griineisen equation of state of the form (@[6l) 

,„ - 73 P - - ^ » ^ I 1 -™ -W ~ ^ - 2V3 r . ; ^:=detF« (4, 

[1 - 6 Q (1 - J C )f 

where Co is the bulk speed of sound, po is the initial mass density, 2 To is the Griineisen's gamma at the 
reference state, S a = dU s /dU p is a linear Hugoniot slope coefficient, U s is the shock wave velocity, U p is 
the particle velocity, and e is the internal energy density (per unit reference volume), F e is the elastic part 
of the deformation gradient. For isochoric plasticity, 

J e = J = det(F) = — . 

P 

The change in internal energy is computed using 

e = PoJ C v dT « p [C V (T) T - C V (T ) T ] 

where To is the reference temperature and C v is the specific heat at constant volume. We assume that C p 
and C v are equal. 

The rate equation (|2]) has three rate terms. We need constitutive relations for each of these three terms. 
The rate form of the pressure equation is obtained by taking the material time derivative of Q to get 

R . Po cl J e [1 + (Sg - 2 r )(l - J e )] ( jA /- 

°™ = ^P= [1 - Sa (l - J*)]3 {Te)- 2V * r ° e 

Now, 

-! = tr (d e ) = 9 e = d e : 1 = V3 d e : 1 ; e = p C V (T) f . 
J e 

For isochoric plasticity, the above relations become 

j r- 

- = tr (d e ) =tr(d) = e = d: 1 = VSd: 1 . 

Hence, 

d m = V3p = 3 K(J e ) d e :1 - 2V3p C v T f 

where 

p C 2 J e [l + (5«-2r )(l- J e )} 



K(J e 



[1 _ Sa (l - J«)]3 

The stress rate a s is given by 



& s = — (VsTs) = s : s = 2 p(a m ,T) rf : s . 

Therefore, 

ci s = 2 p{a m ,T) rf :s . 

The rate of s is obtained by noting that 

s j s : s\ 

s= u\-\m r- 
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Therefore, 

a = [rf - (rj : a) a] . 

A slightly different, but equally valid, decomposition of the stress rate can be obtained by setting 

& s = 2 n(a m ,T) d e :a . 

and 

x 2 fi(a m ,T) r ^ ^ ^; 

s = d — (d : 1)1 — (d : a)a . 

a s l J 

We use the above decomposition in the following development. 



1.3.9 Plastic Relations: 

We consider yield functions of the form 

f(a m , a s , e e p , T) < 

where e is the total strain rate, e p is the plastic strain, and T is the temperature. Here the strain rates and the 
plastic strain are defined as 



e := y | d : d ; e p := \j ^dP : dP ; e p := J e p dt . 

As a particular case, we consider a class of isotropic yield functions of the von Mises form[^] 

tt ■ rp\ ^ 2 2/- rp\ ^ 2 { a mi T) 

f{a m , a s , e, e p , T) = - cr s - cr^e, e p , T) ^ (5) 

where a y is the yield stress, \i is the shear modulus, and \iq is a reference shear modulus. 
At yield, we have 

^ - ^(e e P> T ) 2 = • 

Taking a time derivative of the above equation gives us 

/U 2 day 2 I 1 9fj, 
3a s a s -2a y -, — -2a y -, — a m = 0. (6) 

We assume that the plastic part of the rate of deformation can be obtained from a plastic potential in the 
form (associated flow rule) 

dP = A |£ = A f a 

OCT 

The unit outward normal to the yield surface is given by 

Idf 1 df 
n = 7^- = -zfa ; £ = ||^-|| = /<x h n : n = 1 ==> dP = A f n . 



4 This form of the yield function can be used for deformation induced anisotropic plasticity by converting it to the form used by 
Maudlin and Schiferl (|7|). In addition, the effect of porosity can be incorporated using Brannon's approach (|4|). 
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Therefore, 



'2 . 



o 



Note that 



' 3 

da m - da s 

OCT OCT 



Assuming that e, e p , and T remain unchanged during the elastic part of the deformation, we have 

fa = ^ = ^l+^S = f m l+f s 8; £ = V(/ m ) 2 + (fs) 2 ■ 

da oa m oa s 
For the yield function ([5]), we have 

Jm = -2 O^fe, € p , 1 ) g Q ! Js = 3 (7 S . 



Therefore, equation ([6]) can be written as 



M 2 da v 



f s a s — 2 a y '-^ + fm &rn = ■ 



Also, 



This implies that 



dP = A 



2<(e e P! T) 



^0 da m 



1 + 3 A <To s . 



d p = {dP : 1)1 + (d p : s)s ; d p : ? = -2 A o*(e e p , T) T) ; d p : a = 3 A < 



Recall, 



& = a m 1 + a s s + a s s 
a m = 3 K{J e ) d e :1 - 2^3 p C v T f 
a s = 2 n{<r m ,T) d e : a 
x 2 n(a m ,T) 



d e -{d e : 1)1 - (d e : a)a 



Using the decomposition d e = d — <i p , we get 

& m = 3 K(J e ) {d-d p ):l- 2^3 po C„ T T 
cr s = 2fi(a m ,T) {d-dP) : a 
x 2 n(a m ,T) 



d-dP -{d: 1)1 +(d p : 1)1 - (d:s)s + {dP :s)s 



or. 



& = a m 1 + a s s + a s s 



a m = 3 K(J e ) 



d : 1 + 2Acj2(e, e p , T) 



Mo 



:i <9av, 



2^/3 po C v T t 



& s = 2 p(<7 m , T) (d : s - 3 A a s ) = 2 fi(a m , T) (rj : a - 3 A a s 

a 2fi(a m ,T) 2n(a m ,T) 
s = [rj - (d : a)a] = [»7 - (»7 : *)«] 
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Now, the consistency condition requires that 

m 



f(a m , a s , e e p , T)] = 



or, 
or, 

where 



df &m , df ^ | df de | df . | df f = Q 
<9<r m m <9<r s s <9e <9i de p dT 



fm 6 m + f s & s + / d e + f p e + / t T = 



0/ 5(J v /2d:d ■ „ T l r „ .„ , Tl 



9e p ~ de p ' Jt '~ 8T~ ° v & dT 



Jp - de p ~^ y nl de p ' - «T ~ " ^ iig 

Note that the above equation is the same as equation ([7]). Also note that, 

& m = 3K[d:l-\f m ]- G(T) t ; & s = 2 fi [ V : s - A f s ] 

where 

G(T) := 2V3p C v (T) T . 
Plugging into the consistency condition, we get 

3 K f m [d : 1 - A f m ] - f m G(T) f + 2 /x f s [77 : a - A /,] + f d e + f p e + / t T = 

or, 

A (3 K f* + 2pf s 2 ) = 3K f m d:l + 2 f if a r r .a + f d e + f p e + [f t -f m G(T)} f 

or, 

; _ 3 K f m d : 1 + 2 fi f s -q : s + fd ii + f p e + [ft — f m G(T)} f 

3K/ m 2 + 2/i/ s 2 

1.4 Mass Balance: 

The mass balance equation is 

p + p V • v = 

where p is the current mass density, and p is the material time derivative of p. 

1.5 Momentum Balance: 

The balance of linear momentum can be written as 

V-cr + pb = pv 

where cr is the Cauchy stress, b is the body force per unit volume, and v is the acceleration. 
The momentum equation can be written as 

V • (p 1 +s)+pb = pir 

or, 

V-(pl) + V- s + pb = pv. 
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1.6 Energy Balance: 

The balance of energy can be written as 

a: d— 'V-q + ph = pe 

where q is the heat flux per unit area, h is the heat source per unit mass, and e is the specific internal energy. 
A dot over a quantity represents the material time derivative of that quantity. 

We convert the energy equation into an approximate equation for heat conduction that includes heat source 
term due to plastic dissipation ([SI): 

p C v f - V • (re ■ VT) -(<r:d p = ph 

where C v is the specific heat at constant volume, k is a second order tensor of thermal conductivity, and £ is 
the Taylor-Quinney coefficient. Expressing the stress and the plastic part of rate of deformation into 
volumetric and deviatoric parts, we get 

p C v T- V • (k • VT) - C (p 1 + s) : ( \ tr (d p ) l+rA=ph 
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or. 



P C v T-V-(k- VT) - C Qptr(<F)l : 1 +^tr(d p ) s : 1 + p 1 : rf + s : rf \ =ph 

or, 

p C v f - V • (k ■ VT) - C (p tr (tF) + s : r? p ) = p /i 

2 Weak Forms 

Let us now derive the weak forms of the governing equations. 

2.1 Rate of deformation: 

The weak form of the kinematic relation is given by 

J W* : id - - [Vv + (Vv) T ] | dfl = 

where W* is an arbitrary second-order tensor valued weighting function. We write the rate of deformation 
and the weighting function in terms of their volumetric and deviatoric components to get 

J jdev(W*) + |tr(W*) 1^ : S^q + ^ 9 1 - ^ [Vv + (Vv) T ] j dft = 

or 

^ jdev(W*) de\(W) : 1 - ^de\(W) : [vv + (Vv) 



g -tr(W*) i : r?+ - 0tr(W*) i : i --tr(W*)l : [Vv + (Vv) J J| dfi = 
y |dev(VK*) : J7- ^dev(W*) : [vv + (Vv) T ] + I fltr(W^) - ^ tr(W*) [tr (Vv) + tr ((Vv) T )] | dD. = 
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J jdev(W*) : {r)-d) + ^0tr(W*) - * tr(W*) V ■ v| dfi = 

' "dev(W) : 1 + ^ 0tr(W*) - * tr(W*) V • v 1 dft = 
or 

y | ^ tr (W*) - ^ tr (W*) V • v j dft = 
Define u>* := 1/3 tr (W*). Then the weak form becomes 

y w* (e - v • v) do = o 

2.2 Constitutive Relations: 

The weak forms of the constitutive relations are given by 

y VT** : (a - C(p,T) : r? e ) dft = 

/' 

Jn 



and 

■ P0 c 2 (i-j)[i-r (i- j)] ■ 
p [i - sm - J)} 2 2r ° e 



dn = o 

where W** is a tensor valued weighting function and w** is a scalar weighting function. 



2.3 Momentum Balance: 

To derive the weak form of the momentum equation, we multiply the momentum equation with a 
vector-valued weighting function (w) and integrate over the domain (Q). The weighting function (w) 
satisfies velocity boundary conditions on the parts of the boundary where velocities are prescribed. Then we 
get, 

/ w ■ [V ■ (p 1 ) + V • s + p b] dtt = / pw-irdQ. 
Jn Jn 

Using the identity v • (V • S) = V • (S T ■ v) — S : Vv, where S is a second-order tensor valued field and 
v is a vector valued field, we get (using the symmetry of the stress tensor) 

/ {V ■ [(p 1 ) • w] - p 1 : Vw + V ■ (s ■ w) - s : Vw + p w • b} dQ = / pw-vffl 
Jn Jn 

or, 

/ {V • (p w) — p V • w + V • (s ■ w) — s : Vw + p w • b} dQ = / p w ■ v d£l . 

Jn Jn 

Applying the divergence theorem, we have 

/ {— p V • w — s : Vw + p w • b} dQ + / n ■ (p w) dY + / n ■ (s ■ w) dT = / p w • v dfl 
Jn Jt Jt Jn 

or, 

/ {pw-v + pV-w + s: Vw} d£l — I p w ■ b dQ, + / pn-w dF + / n ■ (s ■ w) dr 

Vn Jn Jt Jt 

or, 

/ {pw-v + pV-w + s: Vw} dfi = / p w • b dfi + / p n ■ w dY + / (s T ■ n) • w dT 

Jn Jn Jt Jt 



15 



/ {pw-v + pV-w + s: Vw} df2 = / p w ■ b df2 + / p n ■ w dr + / (cr T ■ n — p n) • w dr 
in in ir ir 

or, 

/ {pw-v+pV-w + s: Vw} dfi = / pw -bdn+ / (er r • n) • w dT . 
Jn Jn Jr 

In the above, n is the outward normal to the surface T. 

If the applied surface traction is t = cr T • n, we get (since w is zero on the part of the boundary where 
velocities are specified) 

/ {pw-v + pV-w + s: Vw} dQ = p w • b dn + / t • w dT . 

Jn Jn Jr t 

2.4 Energy Balance: 

The weak form of the heat conduction equation is given by 

/ w \p C v f - V • O • VT) - C tr (d p ) + s : r) p )\ dn = [ w phdQ.. 
Jn L i Jn 

Using the identity 0(V ■ v) = V ■ ((f) v) — (V</>) • v we get 

/ w pC v T<m- I { V ■ [w(k ■ VT)] - (Vw) ■ k ■ VT} dQ - f w C [p tr («F) + s : tj"] d!^ = / w p h dQ . 
Jn Jn Jn Jn 

Using the divergence theorem, we have 

J {wpC v T+ (Vw) ■ k ■ VT j dQ, = J w h + C [p tr (d p ) + s : r/ p ] } dfi + ^ n ■ [w?(k ■ VT)] dr . 
If a heat flux (q = n • (re ■ VT)) is specified on part of the boundary (Tg), the weak form becomes 

J \w pC v t + (Vw) ■ k ■ VT} dfl = J w |p h + C [p tr (dP) + s : r? p ]} dn + y to g <fl? . 

3 MPM Discretization 

We can now discretize each of the weak forms using appropriate basis functions. 

3.1 Rate of Deformation: 

The weak form of the rate of deformation equation is 



n 



w* (0 - V • vj dn = . 

We assume a discontinuous approximation for 9 over each cell and a weighting function w* of the form[^] 



p 



#( X ) ~ ^QqXqt*) ! ^*(x) = ^W*Xp( X ) 

<J=1 p=l 



5 The approximations for w* and need not be a sum over the particles. We could alternatively assume a constant value over 
each cell or a sum over the grid nodes with the same basis functions as the grid basis functions. We have chosen a sum over particles 
in a cell to keep our formulation consistent with the standard MPM procedure. 
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where n c p is the number of particles in a grid cell, w* are the weighting functions at the particle locations, 
and Xp( x ) are tne basis functions. The velocities are approximated in the standard MPM way using 



v(x) «£Wx) 

9=1 

where n g is the number of grid points, v 9 are the velocities at the grid vertices, and N g (x) are the grid basis 
functions. 

Plugging the approximations into the weak form, we get the following equations that are valid over each cell 



E™pXp(*) 
P =i 



q=l 



9=1 



dfl = 



{ n P n g 1 

9=1 9 = 1 



The abitrariness of Wp gives us a system of rip equations 
f I" 1 

/ Xp( x ) < 5Z^X 9 (x) - E v s • 

9=1 9=1 



ViV 9 U) = 0; p = 1 ... rip 



£ 

9=1 



Xp(x)X9( x ) 



9=1 



Xp(x)VAT 9 dQ. 



v s = ; p = 1 . . . n v 



9=1 



X P ( x )Xg(x) dfi 



9=1 



X P (x)ViV 9 



v„ ; p = 1 . . . n%. 



If we identify the basis functions Xp with the particle characteristic functions and use the property (J9l, p. 
485) 

I 1 if x £ f2„ 



otherwise 



we get 



/ Xp(x)Xp(x) 



9=1 



v 9 ; p = 1 • • • < 



Define the volume of the particle inside the cell as 

V pc := / x P (x)Xp(x) d£l . 
JflpllQc 

We then have 



Vp C # P - £ 

9=1 

Define the gradient weighting function as 



xp(x.)VN g dn 



■ v g ; p = 1 . . . n 



VN pg :-- 



x P (x)viv 9 dn . 



(8) 
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Then, 



or, 



9=1 



'p - 7^Y, VN P9-vg ; p = i... 

"I - 



9=1 



3.2 Constitutive Relations: 

The weak form of the constitutive relation for pressure is 

p c 2 (i- j(x))[i-r (i- J(x))] 



w**(x) 



p(x) - 



[1 - 5 a (l - J(X))]2 



- 2 r e(x) 



dft = . 



We use a discontinuous approximation for the pressure field (p) and a weighting function w** of the form 



~ ^PgX 9 (x) ; w**(x) = ^u>;*x P (x) . 

g=l p=l 

After plugging these approximations into the weak form, we get the following relations over each cell 



/ 

J si, 



V flv(Tl P oCg(i- J(x))[i-r (i- j(x))] 9Fp , v J rf0 _ n 

l,P^,(x) [l-5 a (l- J(x))P 2roeW dn -° 



P oC,?(i- j(x))[i-r (i- J(x))] 

J(x))]2 



2 r e(x) ) dO, 



Invoking the arbitrariness of w" we get 

f < Jv^ ^ poCgq- j(x))[i-r„(i- j(x))] , r n ) , n n , 

J Q Xp(x) < 2^PlXq(x) [l- g a (l- J(X))] 2 ~ ^ [ = ' P = ' ' 



or, 



Xp(x) X 9 (x) rffi 



/ XpW 



P oC 2 (i- j(x))[i-r (i- J(x))] 

[l-5 a (l- J(x))P 



2 r e(x) 



rffJ ; p = 1 ... rip 



Once again, if we identify the basis functions Xp with the particle characteristic functions and use the 
property 

, 1 1 if x e fi„ 

Xp(x) = 



otherwise 



and we get 



/ XpW Xp(x) dfi 



V'pcPp = / Xp(x) 
Jn„nfi c 



/ 



XpW 



poC 2 (i- j(x))[i-r (i- J(x))] 



[1-S a (l-J(x))]2 



- 2 r e(x) 



p = 1 ... Tip 



PoCgq- j(x))[i-r (i- j(x))] 
[i-s a (i- J(x))] 2 



- 2 r e(x) 



rfri ; p = 1 ... rip 
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Identifying the volume averaged J(x) and e(x) as 9 P and e p , respectively, we get 

p c 2 (i-^)[i-r (i-^)] 



Vpc Pp — Vpc 



[1 - S a (l ~ 8 P )f 



or, 



P P 



P oC 2 (i-^)[i-r (i-^)] 



2 T e p ; p = 1 . . . n c p 



- 2 T e p ; p = \...n c . 



[i-s a (i-e p )Y u p ' 1 •■• p 

The weak form of the constitutive relation for the deviatoric stress rate is given by 



/ 

Jn 



S (x)-C(p,T,x):r, e (x) 



dfi = . 



We assume that the weighting function is 



r*(x) = ^w;* x ,(x) 

9=1 

and plug these into the weak form to get the following relations 

"s(x)-C(p,T,x):*f(x) 



/ 

Jn 



9=1 



dft = 



or, 



£ w 9 * : 

9=1 







/ Xg (x) [l(x)-C(p,T,x):r7 e (x)" 
The arbitrariness of W** gives 

X q (x) s(x) - C(p,T,x) : T7 e (x) d^ = 0; q = l...n p 



J 

Jn 



n a nn 



Using arguments similar to those used for the pressure constitutive relation, we get 



s q = C q (p,T) : rf q ; 

3.3 Momentum Balance: 

The weak form of the momentum balance relation is 



1 . . . n r , 



I [p(x) w(x) • v(x) +p(x) V • w + s(x) : Vw] dQ = p(x) w(x) • b(x) dtt + / t(x) • w(: 
Jn Jn Jr t 

The first step in the MPM discretization is to convert the integral over into a sum of integrals over 
particles. To achieve this we assume that the particle characteristic functions have the form 



Xp(x) = 



1 if x 6 Q p 
otherwise 
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Then the weak form of the momentum balance equation can be written as 

ii p , 

/ [p(x) w(x) • v(x) + p(x) V • w + s(x) : Vw] = 

V / xp(x) p(x) w(x) • b(x) dfi + / t(x) • w(x) dr . 
p=1 ifi p nn ir t 



The weighting function is 



The velocity is approximated as 



w ( x ) = ^2^gN g (x) . 

3=1 



n g 



v(x) ai yjVfejVfe(x) . 

h=l 

The material time derivative of v is also approximated in a similar manner (see |[IOlO . Thus, 

n g 

v(x) « ^v^(x) . 



/t=l 



Plugging these into the left hand side of the momentum equation we get 



LHS = E / Xp( 



p=l " "p 



9=1 



p(x) J ^w 9 iV 9 (x) ■ i £ v fc JV h (x) + p(x) V ■ \ ^w 9 iV 9 (x 



9 = 1 



+s(x): V^w 9 iV 9 (x) 

,9 = 1 



dQ 



LHS = ^\ 



9 = 1 



n -n ( n 



E i E f / Xp (x) p(x) iV 9 (x) iV h (x) dfi ] Vfc 



+ 



E< 

9=1 



J2 XpWp(x) ViV 9 dfl 



9 = 1 



X)/ Xp(x)s(x)- ViV 9 dfi 

p=1 v n p nn 



LHS = ^\ 



9 = 1 



E 



E f / Xp (x) P(x) iV 9 (x) 7V h (x) (01 j v fc 1 + 
h =i \Jn p nn J J 



/ 

Jsn, 



X P (x) p(x) ViV 9 dfi + / Xp(x) s(x) • ViV 9 dfi 



Similarly, the right hand side of the weak form of momentum balance can be written as 



R HS = E / Xp(x)p(x) 



Ew 9 iV 9 (x) 



.9 = 1 



b(x)dfi + / t(x) 



Ls=i 



RHS = ^]w 9 - 

9 = 1 

Tig 

RHS = 1 

9 = 1 



Xp(x)p(x)7V 9 (x)b(x)dfl 

p=1 Jn p nn 



E 1 

9=1 



t(x) JV g (x) dr 



r, 



53 / Xp(x) p(x) JV,(x) b(x) dfi 1 + / t(x) JV,(x) dT 
p=1 JQ p nn J Jr t 
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V ( / Xv (x) p(x) /V 9 (x) N h (x) dfi I v h \ + 
h=1 \Jn p nf! / I 



Combining the left and right hand sides and invoking the arbitrariness of w 9 , we get 



E 

P =i 



V f Xp(x) p(x) 7V fl (x) b(x) dfi I + / t(x) 7V fl (x) dT ; g = 1 . . . n g 

If p(x), p(x), s(x), and b(x) are assumed to be constant over each particle, we can replace these quantities 
with the values at the particles to get 



X P (x) p(x) ViV 3 dft + / x P (x) s(x) • ViV 9 dft 



p=i 



V ( p p ! Xp (x) iV 9 (x) iV ft (x) dQ I v J + 

U=i V ■ /n P nn / J 

, / xp(x) ViV 9 dtt + s p • / x P (x) V7V 9 dO 

Y\ Pp b P / X P (x) iV 9 (x) dft 1 + / t(x) iV fl (x) dT ; = 1... 

Then, using the definition of the gradient weighting function ([8]), and defining 

X P (x) 7V fl (x) dn 



(9) 



we get 







' n p 




' «p "1 


)! 




Y,pp VN p 9 


H 















|jp p b p J + jf t(x) iV 9 (x) dr ; 5 = 1... 



Define the mass matrix as 



m 9h := ^2pp / Xp(x) iVg(x) iVft(x) dfi 
j Jn p r\n 



Define the internal force vector as 



Define the body force vector as 



P =i 



p=l 
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Also define the external force vector as 



Then we get the semidiscrete system of equations 



ff l := / t(x)iV 9 (x)dr 



X> ffh v h + lf = f 9 bod + ff ; 5 = 1. ..n 9 



h=i 



3.4 Energy Balance: 

The weak form of the energy balance equation is 



J {w(x) p(x) C„(x) T(x) + • k(x) • Vr} dft = 

y w(x) {p(x) A(x) + C bW tr (d p (x)) + s(x) : rf (x)]} dft + J w{x) g(x) dT . 
As before, we express this equation as a sum over particles to get 

/ XpW |w(x) p(x) c v (x) t(x) + Vw • «(x) ■ vr} dfi = 

/ Xp(x) u»(x) |p(x) h(x) + C b(x) tr(<F(x)) + S (x) : rf (x)]} dfi + / ffi(x) «(x) dr 
We approximate T using 

n 9 

r(x)«X;r g 7V g (x). 

9=1 

The material time derivative of T is then given by 

r(x)«X;r s iv s (x). 

9=1 

The weighting functions are also chosen to be of the same form: 

^( x ) = ^2^9 N g( x ) ■ 

9=1 

Plugging these into the weak form, we get 

E / X P (x) E^ -°( x ) C -( x ) E^ ^"(x) dQ 

p=1 Jn p nn \ g=1 ) \ h=1 ) 

n p / n g \ / n g \ 

+ E / v E™« ■ «( x ) ■ v E Th ^( x ) dn 

p=1 Jn p nn \ g=1 J \ h=1 ) 

= E / Xp(x) E™« U x ) ^( x ) + C b(x) tr(d"(x)) + S (x) : rf (x)]} dfi 

p=1 ./n p nn \ g=1 ) 1 > 

+ f r (E^^( x )) ^(x)dr 
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For simplicity, let us assume that k(x) is isotropic, i.e., k(x) = «(x), Then, 



E \ E / X P (x) iV 9 (x) p(x) C„ (x) ( ^ iV\ (x) ) dfi 

9 =i Lp=i- /f! f nn u=i / 

+E^iE/ x P (x) «(x) viv a • [Y^tkVNh )<m 

9=1 U^l^p^ \a=1 / 



E- 

9 = 1 



£ / X P (x)7V g (x) {p(x)/ i (x) + C[p(x)tr(d p (x)) + S (x):r ? p (x)]} dtt 



+E*M / ^ 9 (x)g(x)dr 

9=1 lA* 

Invoking the arbitrariness of w g , we get 



]T / Xp (x) iV 9 (x) p(x) C„(x) ( ]T T h JV h (x) ) dn + J2 [ Xp(x) k(x) ViV 9 ■ ( £ T h V7V h ] dfi 
p=1 in p nn \ h=1 / p=1 Ja p nn \ h=1 / 

rip 

= E / Xp(x) JV 9 (x) {p(x) fc(x) + C b(x) tr(rf"(x)) + S (x) : (x)]} dQ 
+ I N g (x)q{x)dF; g = l...n a 



E 



/ xp(x) n 3 (x) iv h (x) P (x) a(x) dn f h +j2 E / ^( x ) *(*) VjV f ■ VNh dQ 



n 



/ Xp(x) N g (x) |p(x) fc(x) + C b(x) tr(d"(x)) + S (x) : r, p (x)]} 



dQ 



+ / iV 9 (x)<?(x) dr ; = l...ra 9 
If we assume that p(x), C„(x), ft(x), h(x), p(x), s(x), and d p (x) are constant over each particle, we get 



E 



/ j Pp Cvp I 



Xp (x)iV a (x)iV fc (x) dfi 



' rip 

I Xp(x) ViV 9 ■ VN h dQ 
p=:1 Jn p nn 



f r 

+ Y,( s p--Vp / X P (x) N g (x) dQ + ATg(x) g(x) dr ; <? = 1 . . . n 9 

p=1 Jn p nfi Jr q 

Using the definition of the particle weighting functions l|9j, we have 



E 



h = i Lp=i 



/] Pp Cv 



X„(x) JV s (x) JV h (x) dfi 



, P =i 



^ k p / Xp(x) ViV 9 • ViVh dfi 



r,, 



^]p P /i P iVp 9 + ^CPptr«) 7V- P 9 + ^ C s P : < iV P9 + / N„{x) q(x) dF ; g = l. 

p=l p=l p=l r 9 



or, 



El Pp / 

p=1 Jn p nn 



X P (x) iV 5 (x) iV fc (x) 



^=1 



E «p / x P (x) viv 9 • VN h an 

j Jn p r\n 



np r i /" 

E K + C P P tr (c^) + C s p : r/Pj Af pg + / 



Af p9 + / iV s (x) ?(x) ; 3 = 1 . . . n g . 
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Define the thermal inertia matrix as 



v 

the thermal stiffness matrix as 



n p 

Mg h :=Y,Pp C vp X P (x) N g (x) N h (x) dfi , 

p=1 Jn p nn 



: =^2 k p ^p( x ) VN 9 ■ VN h dn 

p=1 Jf2 p nn 



the thermal source vector as 

rip 

P =i 

and the external heat flux vector as 



qf x := f N g (x)q( X )dT. 
Then we get the following semidiscrete system of equations: 

ng ri g 

J2 M 9h Th + Y, K 9hT h = S g + qf- g = l. 



n 



9 ■ 



h=l h=\ 



4 Conclusion 



The derivation described in this paper clearly shows that the integration over the volume can be a source of 
errors in MPM. Improved algorithms are required. The stress update algorithm can also be improved by 
using a Lie derivative formulation, possibly with the incorporation of a multiplicative decomposition of the 
deformation gradient. These issues will be addressed in future work. 
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A Useful results 



1 . The integral 



Mt) 

7(t)= f(x,t)dx 

Ja(t) 



is a function of the parameter t. Show that the derivative of F is given by 

,i i /•'"'•' . \ r b(t) df(x,t) 



dt dt 



( fHt) \ fb(t) 

/ /(*,*) dX = / 

\Ja(t) J J a(t) 



dt 



di 



This relation is also known as the Leibnitz rule. 
The following proof is taken from II 11 . 
We have, 



dF F(t + At) - F(t) 

= lim . 

dt At->0 At 



Now, 



F(t + At) - F(t) _ 1 
At ~ At 

1 

= At 
1 

~ At 
1 

~ At 



/ f(x,t + At)dx- / f(x,t), 

Ja(t+At) J a(t) 

/b+Ab rb 
f(x,t + At) dx- / f(x,t)dx 
^+Aa J a 

™+A« pb+Ab rb 

-/ f(x,t + At)dx+ f(x,t + Ai)dK- f(x,t) 

J a J a J a 

/■a + Aa rb rb-\-A.b rb 

- f(x,t + At)dx+ f(x,t + At)dx+ f(x,t + At)dx- f(x,t)dx 

J a J a J b J a 



I 

J a 



b f(x,t + At) - f(x,t) 1 f b + Ab 1 f«+ Aa 

' dx+ — f(x, t + At) dx / f(x, t + At) dx . 

Jb At J a 



At 



At. 
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Since f(x, t) is essentially constant over the infinitesimal intervals a < x < a + An and b < x < b + Ab, we may write 

f 



F(t + At) - FW « f /(M + At)-/(M) dx + /( 6 , t + At) — - /(a, t + At) — . 
At J„ At JK ' ' At Jy ' ' At 



Taking the limit as At — > 0, we get 



lim 

Ai->0 



F(t + At)-F(f) 
At 



lim 

At->0 



J a 



f(x,t + At)~f(x,t) 
At 



dx 



- lim 

Ai->0 



/(M + At) 



Ab 
At 



lim 

Ai->0 



/(a,t + At) 



A a 
At 



at J a(t) at at at 



2. Let Q(t) be a region in Euclidean space with boundary <9f2(f ). Let x(f) be the positions of points in the region and let v(x, t) be the 
velocity field in the region. Let n(x, t) be the outward unit normal to the boundary. Let f (x, i) be a vector field in the region (it may also 
be a scalar field). Show that 



(/ fdv W 

\Jn(t) J Jn 



%dV+ [ (vn)fdA. 

fi(t) Ot JdSl(t) 



This relation is also known as the Reynold's Transport Theorem and is a generalization of the Leibnitz rule. 
This proof is taken from [12] (also see [13]). 

Let Qo be reference configuration of the region f2(f). Let the motion and the deformation gradient be given by 

x = ¥>(X,f); F(X, f ) = Vo . 

Let J(X, f) = det[J r (x, t)]. Then, integrals in the current and the reference configurations are related by 

f f(x,t)dV= / t[<p(X,t),t] J(X,t)dV= / f(X,t) J(X,f) dV . 
in(t) Jn a Jsi 

The time derivative of an integral over a volume is defined as 
f f(x,t), 



d 
dt 



I dV 



lim 

At 



■— - ( f f (x, t + At) dV - / f (x, t) dV ) 

o At yJn( t +At) Jn(t) J 



Converting into integrals over the reference configuration, we get 



d 
dt 



f f(x,t) 



dV 



lim 

At 



oAt(/ ^ X '* + A ') J ( x >* + At)dV-^ f(X,t) J(X,t) dV j . 



Since Qq is independent of time, we have 



,iU (t) f(x < i)dV 



lim 

At->0 



f (X, t + At) J(X, t + At) - f (X, t) J(X, t) 



At 



dV 



n dt 



[f(X,f) J(X,t)] dV 



= jf ^[f(X,t)] J(X,t) + f(X,t) ^[J(X,t)]j dV 

Now, the time derivative of det F is given by (see 1131 , p. 77) 
9J(X,t) _ d 



Therefore, 



dt 
d 

dt 



dt 



(dctF) = (detF)(V • v) = J(X,t) V ■ v(ip(X,t),t) = J(X, t) V • v(x,t) . 



dV 



f f(x,t)dV) = /' f^[f(X,t)] J(X,t) + f(X,t) J(X,t)V-v(x,t) 
Jn(t) J Jn \ot 

= J o (£[f(X,t)] + f(X,t) V-v(x,t)) J(X,t)dV 

= / (f(x,t) + f(x,t) V-v(x,t)) dV 

Jfi(t) v ' 
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where f is the material time derivative of f . Now, the material derivative is given by 

f(x,i) = ^ + [Vf(x, i )].v(x,t). 



at 



Therefore, 



d 

dt \jn{t) 



f f(x, t )dvW (*™ 
Jn(t) I Jn(t) V dt 



dt 



Jn(t) I Jn 



[Vf (x, t)] • v(x, t) + f (x, t) V • v(x, t) ) dV 

) 



( — + Vf ■ v + f V ■ v ) dV 

in(t) \ot 



Using the identity 
we then have 



V • (v (g> w) = v(V • w) + Vv • w 



\Jn(t) J Jnt 



I nit) \9t 

Using the divergence theorem and the identity (a ® b) • n = (b • n)a we have 
d I f „ . \ f df 



V • (f <g> v) dV 



f f dV ) = f ^ dV + f (f ® v) • n dV = / 

Jn(t) I Jn(t) dt Jdn(t) Jn(t) 



dt \Jn{t) 

3. Show that the balance of mass can be expressed as: 



df 

di 



dV 



an(t) 



(v • n)f dV . 



p + p V • v = 



where p(x, t) is the current mass density, p is the material time derivative of p, and v(x, t) is the velocity of physical particles in the body 
SI bounded by the surface dil. 



Recall that the general equation for the balance of a physical quantity /(x, t) is given by 

/ /(x,t)dV =/ /(x,t)[u n (x,t)-v(x,t)-n(x,t)]dA+ / g(x,t)dA+/ fc(x,t)dV. 
Jn J Jan Jan Jn 



d 
dt 



To derive the equation for the balance of mass, we assume that the physical quantity of interest is the mass density p(x, i). Since mass is 
neither created or destroyed, the surface and interior sources are zero, i.e., g(x, t) = /i(x, t) = 0. Therefore, we have 



d 
dt 



f p(x,t)dV = / p(x,t)[u„(x,t) - v(x,t) • n(x,t)] dA. 
Jn J Jan 



Let us assume that the volume SI is a control volume (i.e., it does not change with time). Then the surface dfl has a zero velocity («„ = 0) 
and we get 

I % AW = - I P ( v ' n ) ' 

Jn at J gn 



I dA . 



Using the divergence theorem 



we get 



v • n dA 



>an 



f V vdV= / 

Jn J a 

V-(pv) 



dV. 



f \°p 

Jn 



Since S7 is arbitrary, we must have 

Using the identity 
we have 

Now, the material time derivative of p is defined as 
Therefore, 



dt 

dp 
at 



dV = 0. 



- V • (pv) = . 



V • (<p v) = ip V • v + V<p • v 

dp 



dt 



+ p V • v + Vp • v = . 



p = — + Vp ■ v . 

dt 



p + p V • v = . 



27 



4. Show that the balance of linear momentum can be expressed as: 

p v — V<t — pb = 

where p(x, t) is the mass density, v(x, t) is the velocity, <r(x, t) is the Cauchy stress, and p b is the body force density. 
Recall the general equation for the balance of a physical quantity 



d 
dt 



f /(x,t)dV =/ /(x,t)[u n (x,t) - v(x,t) ■ n(x,t)] dA+ f 9 (x,t)dA+/ ft(x,t)dV. 
in J Jan Van in 



In this case the physical quantity of interest is the momentum density, i.e., /(x, t) = p(x, i) v(x, t). The source of momentum flux at the 
surface is the surface traction, i.e., 9 (x, t) = t. The source of momentum inside the body is the body force, i.e., 
/i(x, t) = p(x, t) b(x, t). Therefore, we have 



d 
dt 



j pvdV = / p v[u n - v • n] dA+ / tdA+ / pbdV 

in J Jan Jan Jn 



The surface tractions are related to the Cauchy stress by 
Therefore, 



t = cr • n . 



d 
dt 



/ p v dV = / p v[u n - v • n] dA + / cr • n dA + / p b dV . 

.Jn J Jan Jen Jn 

Let us assume that Q is an arbitrary fixed control volume. Then, 

f ^(pv)dV = - f pv(v-n)dA+ / <r • n dA + /" p b dV . 

Jn dt Jgn Jan Jn 

Now, from the definition of the tensor product we have (for all vectors a) 

(u ® v) • a = (a • v) u . 



Therefore, 



Using the divergence theorem 



we have 



= - / P (v ® v) • n dA + / cr-ndA+/pbdV. 

Jan Jan Jn 

/ V • v dV = / v ndA 

Jn Jan 

/ ^-(pv)dV = - /" V • [p (v® v)] dV + f V <rdV+ f p b dV 

Jn ot Jn Jn Jn 

In [^ (pv) + V ' [(/9v) ® vl_ V '°--P b 



dV = . 



Since is arbitrary, we have 

Using the identity 

we get 



at 



(p v) + V • [(p v) ® v] - V • cr - p b = . 



V • (u <g> v) = (V • v)u + (Vu) ■ v 



1ft V + P ^ + (V ' V)(pV) + V (' ov )- v - V -' T -' :,b = 



rip 



+ pV-v 



9v 

v + p h V(p v)-v-V-cr-pb = 

at 



Using the identity 

we get 

From the definition 



dp 
dt 



+ pV-v 



V(<p v) = <p Vv + v ® (V<p) 



v + p — + [p Vv + v (g> (Vp)] •v-V-cr-pb = 



(u ® v) • a = (a • v) u 
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we have 



Hence, 



dp 

at 



p V ■ v 



[v <g> (Vp)] • v = [v (Vp)] v . 



dv 



v + p — + p Vv • v + [v • (Vp)] v-V<r-pb = 



dp 
at 



Vp • v + p V ■ v 



dv 

v + p h p Vv v — V ■ <r — pb = 0. 



The material time derivative of p is defined as 



Therefore, 

From the balance of mass, we have 
Therefore, 

The material time derivative of v is defined as 
Hence, 



dp _ 
p = — + Vp • v . 

dt 



9v 



[p + p V • v] v + p — + p Vv v — V • <r — pb = 0. 



p + p V • v = . 



dv 



p h p Vv v — V • <r — pb = 0. 

9t 



<9v 

at 



Vv • v . 



p v — V • <r — pb = 0. 
5. Show that the balance of angular momentum can be expressed as: 



We assume that there are no surface couples on c9S! or body couples in fl. Recall the general balance equation 

f /(x,t)dV =[ /(x,t)[i t „(x,t)-v(x,t)-n(x,t)]dA + / 9 (x,t)dA+/ ft(x,t)dV. 
_Jn J Jan van Jn 



(it 



In this case, the physical quantity to be conserved the angular momentum density, i.e., / = x X (p v). The angular momentum source at 
the surface is then j = xxt and the angular momentum source inside the body is h = x X (p b). The angular momentum and moments 
are calculated with respect to a fixed origin. Hence we have 



d 
dt 



/ x X (p v) dV = / [x X (p v)] [tt„ - v • n] dA + / x X t dA + / x X (p b) dV . 

Jn J Jan Jen Jn 



in J Jen 

Assuming that Q is a control volume, we have 
' d 



Jn 



dt 



(pv) 



dV = - / [x X (p v)] [v • n] dA + / x X t dA + / x X (p b) 

Jan Jan Jn 



dV. 



tan Jan 

Using the definition of a tensor product we can write 

[x X (p v)] [v • n] = [[x X (p v)] <g> v] • n . 

Also, t = <t • n. Therefore we have 
' d 



( 

Jn 



x x 



(PV) 



dV = - f [[x X (p v)] <g> v] • n dA + f x X (<r • n) dA + / x X (p b) dV 



Using the divergence theorem, we get 
' d 



Jn 



at 



(pv) 



dV = - / V • [[x X (p v)] ® v] dV + / x X (<r • n) dA + / x X (p b) 

Jan Jn 



dV . 



To convert the surface integral in the above equation into a volume integral, it is convenient to use index notation. Thus, 

yx X (er • n) dA = / e ijfe Xj o kl n ; dA = / An n\ dA = / A • n dA 
en Ji Jan Jan Jan 
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where [ ]; represents the i-th component of the vector. Using the divergence theorem 



f A-ndA= f V-4dV= f dV = f — (e ijk x 3 - a kl ) dV . 

Jan Jn Jn dxi Jn dxi 



Differentiating, 



f AndA= f 

Jan Jn 



^ijk &jl ^kl ~t~ &ijk X 



9<T k l 

3 dxi 



dV : 



Jn . 



^ijk ^kj &ijk %j 



dx. 



dV 



/ [eijk °kj + e-ijk %j [V • cr]i] dV . 
Jn 



Expressed in direct tensor notation, 



/ A ■ n dA = f \[£ : <r T ] l + [x X (V • cr)]i] dV 
Jan Jn L J 

where £ is the third-order permutation tensor. Therefore, 



xx (cr • n) dA 



[[£ : <r T ] l + [x X (V • dV 



i an Jn 
The balance of angular momentum can then be written as 



/ x X (cr • n) dA == / \s : cr T + x X (V • cr)] dV . 

Jan Jn 1 



f 

Jn 



x x 



(PV) 



dV = - / V • [[x X (p v)] ® v] dV + / [f : <r T + x X (V • <t)1 dV + / x X (p b) dV . 
Jn Jn 1 J Jn 



Since f2 is an arbitrary volume, we have 
' d 



x X 



at 



(pv) 



- V • [[x X (p v)] ® v] + £ : cr T + x X (V • cr) + x X (p b) 



x X 



— (pv) - V • cr - pb 



-V-[[xX (p v)] ® v] + £ : cr 1 



Using the identity, 
we get 



V • (u <g> v) = (V • v)u + (Vu) • v 



V • [[x X (pv)] ® v] = (V • v)[x x (pv)] + (V[x x (pv)]) • v . 
The second term on the right can be further simplified using index notation as follows. 



d 

[(V[x X (pv)])-v], = [(V[p (x X v)])- v] t = —(pe ijk Xj v k ) v t 



' e ijk 



Xj V k Vt + p V k V i + p Xj VI 

oxi 



dxi 



lJk [dxi 3 

= (eijk ^j Vk) + P ( e ijk S jl V k v l) + e ijk Xj ^p V^j 

= [(x X v)(Vp -v) + pvXv + xX (p Vv • v)]» 
= [(x X v)(Vp • v) + x X (p Vv • v)]j . 

Therefore we can write 

V • [[x x (p v)] ® v] = (p V • v)(x X v) + (Vp • v)(x x v) + x x (p Vv • v) . 
The balance of angular momentum then takes the form 



x X 



a 

— (p v) - V • <t - pb 



-(p V • v)(x X v) - (Vp • v)(x X v) - x X (p Vv • v) + £ : cr 1 



x X 



at 



(pv) + p Vv • v — V • cr — p b 



-(p V ■ v)(x X v) - (Vp • v)(x X v) + £ : cr 1 
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x X 



dv dp „ „ 

+ — v + pVvv-Vcr-pb 



-(p V • v)(x X v) - (Vp ■ v)(x X v) + £ : er T 



The material time derivative of v is defined as 



dv 

v=-+Vv.v. 



Therefore, 



xx [p v - V • <r - p b] = -x X ^ v + -(p V • v)(x X v) - (Vp • v)(x X v) + £ : cr T . 



Also, from the conservation of linear momentum 
Hence, 



pv — Vcr — pb = 0. 



= x X ^ v + (p V • v)(x X v) + (Vp • v)(x X v) - £ : cr T 
at 

dp 
di 



pV • v + Vp • v I (x X v) - £ .a- 1 



The material time derivative of p is defined as 



dp _ 
p= _ +V p.v. 



Hence, 

From the balance of mass 

Therefore, 

In index notation, 

Expanding out, we get 

Hence, 



(p + p V • v)(x X v) - £ : cr T = . 
p + p V • v = . 
£ : (T T = . 

£ijk crfcj = . 
""12 — o"2i = ; (723 — C32 = ; (731 — (713 = . 



6. Show that the balance of energy can be expressed as: 

p e — <t : (Vv) + V • q — p s = 

where p(x, t) is the mass density, e(x, t) is the internal energy per unit mass, <r(x, t) is the Cauchy stress, v(x, t) is the particle velocity, 
q is the heat flux vector, and s is the rate at which energy is generated by sources inside the volume (per unit mass). 

Recall the general balance equation 



d 
dt 



f /(x,t)dV = f /(x,t)[tt„(x,t)-v(x,i)-n(x,t)]dA+ / <?(x,t)dA+/ h(x,t)dV. 
in J Jen Jan Jn 



In this case, the physical quantity to be conserved the total energy density which is the sum of the internal energy density and the kinetic 
energy density, i.e., / = p e + 1/2 p |v • v|. The energy source at the surface is a sum of the rate of work done by the applied tractions 
and the rate of heat leaving the volume (per unit area), i.e, g = v ■ t — q ■ n where n is the outward unit normal to the surface. The energy 
source inside the body is the sum of the rate of work done by the body forces and the rate of energy generated by internal sources, i.e., 
h = v • (pb) + p s. 

Hence we have 



d 
dt 



f 



p I e H — v • v ) dV 



f p ( e + i v • v ) (u n - v • n) dA + f (v • t - q • n) dA + f p (v • b + s) dV . 

Jan \ 2 J Jqq Jn 



Let £7 be a control volume that does not change with time. Then we get 



L 



d_ 
a dt 



f j i e+ i v- v 



dV : 



/ p (e+-v-v) (v-n)dA+ f (v • t - q • n) dA + f p (v • b + s) dV . 

Jan V 2 ) Jqq Jn 



/an \ / Jan 

Using the relation t = a ■ n, the identity v • (er • n) = (<r T • v) • n, and invoking the symmetry of the stress tensor, we get 

In a\ \ p ( e + \ V ' V ) dV = _ f g ( e + ^ v ' v ) ( v ' n ) dA + ' v ~ q) ' n dA + ^ -° ( v ' b + s ) dv ■ 
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We now apply the divergence theorem to the surface integrals to get 



/ 











dV = - / V • 


P \e + \ v • v 









dA + / V • (<7 • v) dA - / V • q dA + / p (v • b + s) dV . 



Since f2 is arbitrary, we have 



_9 

at 



p ( e+ - v v 



P \ e + 2 v ' v l v 



+ V • (<r • v) - V • q + p (v • b + s) . 



Expanding out the left hand side, we have 



d_ 
dt 



( 1 M d P ( 1 ^\ (de 1 d . ,\ 
p ^+-vvjj = _^ + - v .vj+p ^_ + -_(v.v)j 



(9/9 



at 



c9e 9v 



— e+-v-v +p — + p — - v 



dt St 



For the first term on the right hand side, we use the identity V • (ip v) = <p V • v + Viy3 ■ v to get 



p ( e + - v ■ v ) V-v + V 



p \ e+ - v v| v 



p ( e+ - vv 



p (e+ — v-v) V-v 



p ( e + - v ■ v ] V-v 



v ■ v ) Vp - v + pV(e+ — v-v) - v 



v ■ v ] Vp - v + pVe-vH — p V(v ■ v) • v 



: p + - v • v^ V • v + + - v ■ v^ Vp • v + p Ve • v + p (Vv 1 



p ( e + - v • v ) V-v 



1 



v-v Vp ■ v + p Ve - v + p (Vv • v) - v . 



For the second term on the right we use the identity V - (S T ■ v) = S : Vv + (V • S) ■ v and the symmetry of the Cauchy stress tensor 
to get 

V - (<r - v) = <T : Vv + (V - <r) • V . 



After collecting terms and rearranging, we get 



'dp 



1 



p V • v + Vp - v) leH — v ■ v I + I p — — h p Vv • v — V • cr — p b 



dv 



dt 



3 ). v + p (| + Ve.v) 



- <T : Vv + V q-ps = 0. 



Applying the balance of mass to the first term and the balance of linear momentum to the second term, and using the material time 
derivative of the internal energy 

de 

we get the final form of the balance of energy: 

p e — cr : Vv + V- q — p s = . 
1. Show that the Clausius-Duhem inequality in integral form: 

(/ pr?dV) > / pr,(n n - vn)dA- / ^dA+Z^dV. 

\Jn J Jan Jan J- Jn J- 



d 
dt 



can be written in differential form as 



+ ■ 



ps 



pr)> -V • 

Assume that ft is an arbitrary fixed control volume. Then u n = and the derivative can be take inside the integral to give 
f ^(pr?)dV>- / pjj(vn)dA- / - 1 — " dA + [ — dV . 
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Using the divergence theorem, we get 

ps 



/«^^ dV ^-I V( ^ V)dV -/« V (l) dV + /o 



dV. 
T 



Since Q is arbitrary, we must have 



^( p7? )>-V-(pr,v)-V. (|) 



Expanding out 



l^g.-v^.v-^v-v)-*. 3 + £i 



dp dr\ / q \ p s 

^■') + P^>-')V f )-v-pVt)-v-pi l (V-v)-V- - +y 



'| + Vp.v + pV.v) , + p (*Z +V „.v)>-V.(|j+^. 



Now, the material time derivatives of p and J7 are given by 

dp d-q 
p = — + Vp ■ v ; ?? = — + V»7 • v . 

at at 

Therefore, 

/ q\ p s 

(p + p V • v) r) + pr)> V • 
From the conservation of mass p + p V • v = 0. Hence, 



T I T 



8. Show that the Clausius-Duhem inequality 

n > —V • I 

, T I T 

can be expressed in terms of the internal energy as 

p{e-Tri)-<r: Vv< — . 

Using the identity V-(ipv) = <pV-v + v- V<p in the Clausius-Duhem inequality, we get 

(q\ ps 1 / l\ p s 

0r M>-yV-q-q-V -J +— . 

Now, using index notation with respect to a Cartesian basis ej , 

A 9 i -i , m 9T 1 „ m 
- = ) e, = - (T ) e,- = VT . 

Hence, 

^>-Jv.q+^q.VT+^ or p ^ > - 1 (V • q - p s) + -L q ■ VT . 
Recall the balance of energy 

p e — cr : Vv + V ■ q — p s = p e — er : Vv = — (V ■ q — p s) . 

Therefore, 

1 , 1 q • VT 

pfl> -(/)e-(r:Vv)+- q-VT => pr]T>pe~cr: Vv H — — . 

Rearranging, 

, q-VT 

p (e - T r)) - a- : Vv < — . 
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9. For thermoelastic materials, the internal energy is a function only of the deformation gradient and the temperature, i.e., e = e(F, T). 
Show that, for thermoelastic materials, the Clausius-Duhem inequality 



p(e-Tf,)-tr: Vv< — 

can be expressed as 

de \ ( de „ T \ . q ■ VT 

T ) n+ p cr - F~ T : F < . 

dr, ) 1 V 9F J ~ T 

Since e = e(F, T), we have 

de ■ de 

e = : F -\ n . 

<9F 9r; ' 

Therefore, 

/ de . de m \ „ q • VT f de \ de . „ q • VT 

'(^♦o^-^)-* 5 ^- — ° r "Ur T ) 7, + p-:F- CT :Vv<-^. 

Now, Vv = I = F ■ F -1 . Therefore, using the identity A : (B ■ C) = (A ■ C T ) : B, we have 

cr : Vv = cr : (F ■ F" 1 ) = (cr • F~ T ) : F . 



Hence, 



, de \ de ■ , „ T , ■ q • VT 

p I — T ) r) + p — — : F — (cr ■ F ) : F < 



dr) J r dF T 

q- VT 



'(£- T M'£-'-'i = - 



10. Show that, for thermoelastic materials, the balance of energy 

p e — cr : Vv + V • q — p s = . 



can be expressed as 
Since e = e(F, T), we have 

Plug into energy equation to get 

Recall, 

Hence, 



pTtj = -V-q + ps. 



de ■ de 
e = : F H ?? . 



de ■ de 
p — : F + p — 7)-cr:Vv + Vq-ps = 0. 
dF dr] 



de m , <9e „ t 

— = T and p = cr • F" T . 

9t; P 9F 



(cr • F~ T ) :F + pT7)-cr:Vv + V- q- ps = 0. 
Now, Vv = I = F ■ F _1 . Therefore, using the identity A : (B • C) = (A ■ C T ) : B, we have 

cr : Vv = cr : (F • F" 1 ) = (cr • F~ T ) : F . 

Plugging into the energy equation, we have 

cr : Vv + p T r) - cr : Vv + V- q- ps = 

or, 

pT r) = —V • q + p s . 
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1 1 . Show that, for thermoelastic materials, the Cauchy stress can be expressed in terms of the Green strain as 



<r = pF ■ 



dE 



Recall that the Cauchy stress is given by 



<r = p • F 1 

F dF 



Oij = P ■ 



de 



dF ik ~ kj P dF lk 

The Green strain E = E(F) = E(U) and e = e(F, r]) = e(U, r]). Hence, using the chain rale, 



de 



de _ de dE 
dF ~ dE ' dF 



Now, 



E = -(F T ■ F - 1) 



de de dE lm 

dF ik dE lm dF ik 

Elm = 7)( F lp F vm ~ <5; m ) = -(F p i F } 



Taking the derivative with respect to F, we get 



dE _ 1 f dF 
dF ~ 2 



V dF dF) 



dEi m 1 / dF p i dF pm 

1 fpm + r v i 



dF ik 2 \dF lk 



dF tl 



Therefore, 



1 



de / dF T 
2 r [dE ' V dF 



■ F + F 



Recall, 



dF\ 
d~Fj\ 

dA = dAij 
dA = dAki 



n 3 = -P 



de fdF pl 



dF„ 



dE lm \dF t 



dFit 



<>ik Sji and 



Therefore, 



de 



dE lm 

de 
dE km 



de . 



dA T _ dAji 
-dA=dA k - l - 5 > k5u 

de 



1 



■jk 



2 VdE ln 



(Sit + F U S mk ) 



"jk 



dE lk 



fjk 



a = - p 

2 P 



F (I) 



T de 
+ F -dE 



p jk 



<r=-pF- 

2 



/ de 



X 



de 



\dE J dE 



From the symmetry of the Cauchy stress, we have 



cr = (F ■ A) ■ F T and cr T = F ■ (F ■ A) T = F ■ A T ■ F T and <r = cr T =$> A = A 1 



Therefore, 



and we get 



de _ / de\ T 
dE ~ \dEJ 

„ de „ 7 

a = pF F 

1 dE 



12. For thermoelastic materials, the specific internal energy is given by 

e = e{E,rj) 

where E is the Green strain and t) is the specific entropy. Show that 

— ( e -Tr)) = -fr)+—S:E and — (e - T r, - — S : E) = -f r) - — S : E 

dt po dt po po 

where po is the initial density, T is the absolute temperature, S is the 2nd Piola-Kirchhoff stress, and a dot over a quantity indicates the 
material time derivative. 
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Taking the material time derivative of the specific internal energy, we get 

de ■ de 

e = : E H r? . 

dE dt) ' 

Now, for thermoelastic materials, 

m de . „ de 

T = — and S = po . 

dt) dE 

1 1 
e=—S:E + Tf). => e-Tf]=—S:E. 
PO po 

^-(Tr 1 )=fri + Tri. 
at 

e-^{Tr,)+tr,= —S:E => ^-{e — T-q) = —T r) + — S : E . 
at po at po 

i \ i . i . 

. —S:E\ = —S : E + —S : E . 

dt \ po I po po 



Therefore, 



Now, 



Therefore, 



Also, 



Hence, 



e-^(Tri)+f V =*-(—S-.E)-—S-.E => ^ I e - T r, - — S : E I = -T V - — S : E . 
dt dt yp J p dt \ po I po 

13. For thermoelastic materials, show that the following relations hold: 

^ = is(£,T); ^ = -fj(E,T) ; ^ = ^E(S,T); ^-=fj(S,T) 
dE po v ' dT ly " dS po 9T ,v ' 

where V>(-^> T) is the Helmholtz free energy and g(S, T) is the Gibbs free energy. 
Also show that 

dS drj dE dfj 

— = - P0 — and ^=P0^- 

Recall that 

iP(E,T) = e-T V = e(E,r l )-T V . 

and 

g(S, T) = -e + T V + —S : E . 

PO 

(Note that we can choose any functional dependence that we like, because the quantities e, r\, E are the actual quantities and not any 
particular functional relations). 

The derivatives are 

dtp _ de _ 1 dtp _ 

d~E~ d~E~ ~p~o~ ' dT ~ ~ T> ' 

and 

dg 1 dS 1 dg 

— = :E=—E; — = n . 

dS po dS po dT ' 

Hence, 

°± = ±S(E,T); ^=-v(E,T); ^ = ^E { S,T); ^=v(S,T) 
dE po dl do po ol 

From the above, we have 

d 2 tp _ d 2 tp _ dfi _ 1 dS 

dTdE ~ dEdT ~aE ~ po<9T ' 

and 

d 2 g _ d 2 g dfj _ 1 dE 

dTdS ~ dSdT ^ dS ~ ~po~&T ' 

Hence, 

dS dfj , dE dfj 

= — on and = po . 

dT F dE dT F dS 
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14. For thermoelastic materials, show that the following relations hold: 



and 

Recall, 
and 

Therefore, 
and 

Also, recall that 



de{E,T) _ df) _ dH 
dT dT dT 2 



de(S,T) _ T df L + ^ s dE =T (Pi + 5 d 2 g 



dT dT po dT dT 2 dSdT 

i,(E,T) = Tp = e-Tr, = e(E,T)-Ti)(E,T) 

g(S, T)=g = -e + T V +—S:E = -e(S, T) + T fj(S, T) + — S : E{S, T) 
Po po 

de(E,T) _ dj, T) , T 9V 
dT -df +v{E ' T) + T df 

de(S,T) 9g di) 1 dE 

dT = _ dT T) +T — + —S . 

V(E,T) = 
V(S,T) 



and 



Hence, 



and 



dip 




df) 


~df 




df 


_ dg 




df) 


dT 




df 


dg 
' dS 




dE 

df 



,)■•,, 



de(E,T) _ df, _ T dH 
dT dT dT 2 



de{S,T) _ T df L + ^ s dE =T & 2 l + s d 2 g 



dT dT po ' 9T dT 2 ' dSdT 

15. For thermoelastic materials, show that the balance of energy equation 

p T fj = — V • q + p s 



can be expressed as either 



where 



p C v f = V • (« • VT) + p s + — T || : E 

Pa dT 



» \ C P --S:^) f = V-(K-VT)+ps-^T^:S 
Po oT I po dT 



C v = d A^l and C p ^ g(5 ' T) 
dT 1 



If the independent variables are E and T, then 



v = f)(E,T) => r, = ^..E + ^f. 
' IK ' ' ' dE dT 



On the other hand, if we consider S and T to be the independent variables 



Since 



v = fj(S,T) => ^^-.S+^f. 

' dS dT 



df) _ _ 1 dS df) _C V dfj _ 1 dE dfj _ 1 / _ 1 dE 

dE ~ ' 'p^df ' df ~ ' ~T " ' dS ~ po df ' df ~ f\ p ~ pf, ' df 
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i as ■ c v . 

f) = : E H T 

po dT T 

1 dE if 1 dE\ ■ 

f)= — — :S+-[C P S:— T . 



we have, either 



PO dT Ty po <9T y 

The equation for balance of energy in terms of the specific entropy is 

pT f) = — V ■ q + p s . 
Using the two forms of f), we get two forms of the energy equation: 

- — :E + pC v f = -V-q + ps 

PO oT 



and 



p BE ■ p dE ■ 

— T :S + pC v T-—S: T = - V ■ q + p 

Po dT po dT 



and 



Rearranging, 



From Fourier's law of heat conduction 

q= —k • VT . 

Therefore, 

-^T^:B + pC„t = V-( K -VT) + f)S 
PO 

p dE ■ p dE ■ 

— T — : S + pC p T - JL S : — T = V • it • VT) + p s . 

PO «T po oT 

p C„ T = V ■ (re • VT) + p s + — T ^ : .E 

Po oT 

' ; I C p~ ~ S: f = V-(K-VT)+p S - -T^ :S. 

po oT y po dT 

16. For thermoelastic materials, show that the specific heats are related by the relation 

Cry — Cv — I S — T I '. . 

P Po V dT J dT 

Recall that 

_ de(E,T) _ T df, 
dT dT 



c _ de(S,T) _ T df, | 1 g dE 



and 

r '" 

<9T dT po dT 

Therefore, 

C P -C„=T^ + 1S:^-T^. 
P dT po dT dT 

Also recall that 

U = »?(E,T)=ij(S,T). 
Therefore, keeping S constant while differentiating, we have 

df\ df) dE df) 
df ~ dE ' df + df ' 

Noting that E = E(S, T), and plugging back into the equation for the difference between the two specific heats, we have 

p dE dT po dT 

Recalling that 

drj _ _ 1 dS 
dE ~ ~pf,df 

we get 

p po \ dT BT 
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17. For thermoelastic materials, show that the specific heats can also be related by the equations 

q2„ 



P V ~po ' df + df ' \8EdE ' df) ~ ~p~o ' ~9~f + Po &T ' \~dE~ ' df ) 

Recall that 

S = po = po f(E(S,T),T) . 

Recall the chain rule which states that if 

g(u,t) = f(x(u,t),y(u,t)) 
then, if we keep u fixed, the partial derivative of g with respect to t is given by 

dg df dx df dy 
dt dx dt dy dt 

In our case, 

u = S, t = T, g(S,T) = S, x(S,T) = E(S,T), y(S,T) = T, and / = p /. 

Hence, we have 

S = g(S, T) = f(E(S, T), T) = po f(E(S, T), T) . 
Taking the derivative with respect to T keeping S constant, we have 

11 





dT fff 



Po 



df dE df &f 
dE ' df + df fff 



Now, 



d^.dE.dl 
dE ' dT dT 



^ dip df d 2 <p df d 2 iP 

f = — => t~= = and 



Therefore, 



Again recall that, 



dE dE dEdE dT dTdE 

d 2 tp dE d 2 ip d (dip\ dE d / dip \ 
~ dEdE ' df + dTdE ~ dE \dE ) ' df + df \dEJ 

9± = ^ s 
dE po 

Plugging into the above, we get 

d 2 ip dE 1 dS _ 1 dS dE 1 dS 
° ~ dEdE ' df ' + ~p7,df ~ po~dE ' df ' + p~o~df ' 

Therefore, we get the following relation for dS /dT: 

dS _ _ d 2 ip 8E _ _dS dE 
df ~ ~ Po dEdE ' df ~ ~ dE ' df ' 



Recall that 

C-n — Oii — — I S — T I ; 

dT I dT 



Plugging in the expressions for dS /dT we get: 



po V 



d 2 ip dE\ dE 1 / dS 8E\ dE 

S + Tpo — ^- :^ :^ = — [S + T 



1 po V dEdE dT ) dT po \ dE dT J dT 

Therefore, 

1 dE / d 2 ip dE\ dE _ 1 dE T / dS 8E\ dE 

p ~ " ~ pf, 'df + \dEdE : df) : df ~ pf, ~df + pf, \d~E ' Iff ) ' Iff 

Using the identity (A : B) : C = C : (A : B), we have 

c ~c =-s- — +t— - (J^L .fw\ = L s -— + — — •(—•—') 

P po dT dT ' \dEdE ' dT J po ' dT po dT ' \dE ' dT ) 
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Consider an isotropic thermoelastic material that has a constant coefficient of thermal expansion and which follows the St-Venant 
Kirchhoff model, i.e, 

a E = a 1 and C = Ai(gll+2/il 
where a is the coefficient of thermal expansion and 3A = 3K — 2/i where K, fi are the bulk and shear moduli, respectively. 
Show that the specific heats related by the equation 

C p - C v = — \a tr (S) + 9 a 2 K T] . 
PO 

Recall that, 

1 T 

C p — C v = — S : a E H ol e ■ C : OL E ■ 

PO PO 

Plugging the expressions of ole an d C into the above equation, we have 

1 T 
C p - C v = — S : (a 1 ) H (a 1 ) : (A 1 ® 1 + 2/i I) : (a 1 ) 

PO PO 

a a 2 T 
= — tr (5) + 1 : (A 1 ® J + 2// I) : 1 

po PO 

a a 2 T 
= — tr (S) + 1 : (A tr (i ) 1 + 2/u i ) 

PO PO 

= — tr(S) + C ^ r (3Atr(J) + 2 M tr(J)) 
Po Po 

a , s 3 a 2 T 

= — tr(S) + (3A + 2 M ) 

Po Po 

atr(S) 9a 2 iv"T 

= —+ . 

Po Po 



Therefore, 



— \atr(S) + 9a 2 KT] 
PO 
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